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HIGH-ORDER “CYCLO-DIFFERENCE” TECHNIQUES: AN ALTERNATIVE TO 

FINITE DIFFERENCES 


Mark H. Carpenter * and John Otto * 


ABSTRACT 


The summation-by-parts energy norm is used to establish a new class of high-order finite-difference 
techniques referred to here as “cyclo- difference” techniques. These techniques are constructed cycli- 
cally from stable subelements, and require no numerical boundary conditions; when coupled with the 
simultaneous approximation term (SAT) boundary treatment, they are time asymptQtically stable for 
an arbitrary hyperbolic system. These techniques are similar to spectral element techniques and are 
ideally suited for parallel implementation, but do not require special collocation points or orthogonal 
basis functions. The principal focus of this work is on methods of sixth-order formal accuracy or less; 
however, these methods could be extended in principle to any arbitrary order of accuracy. 

INTRODUCTION 


A great deal of effort has recently been placed on high-order finite-difference techniques (both 
central and upwind) for direct numerical simulations. A significant problem that faces the high-order 
finite-difference community is the closure of those boundary schemes that retain the formal accuracy 
of the underlying method and do not cause instability. To retain the Nth - order formal accuracy of 
the interior scheme for an arbitrary hyperbolic equation, the numerical boundaries must be closed 
with an accuracy of no less than ( N — 1 )th order [1]. Such closures often cause numerical instability 
and cannot be used, (e.g., [2]). Recently, a precise means of determining boundary closures that 
maintain both stability and accuracy has been developed based on the summation-by-parts energy 
norm. (See references [3], [4], and [5].) A numerical discretization that satisfies specific criteria on the 
discretization matrix A * automatically satisfies a discrete energy norm. Central-difference schemes 
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automatically satisfy these properties in their interior. The task is, therefore, to find high-order 
closure techniques at the boundaries that maintain the specific form of the matrix A*. This can be 
a daunting task, but can be accomplished [2], [4], [5]. 

One of the advantages of spectral techniques is that the ambiguity of numerical boundary treat- 
ments is not present. The global nature of the method relies on a specific stencil at each point. Each 
individual point may be unstable, but the scheme as a whole is stable and accurate. This notion 
motivates us to develop a new class of finite-difference schemes. Like central-difference and spectral 
techniques, they are not biased in the direction of a physical eigenvalue and, therefore, do not re- 
quire eigenvalue decomposition when used for general hyperbolic equations. Unlike central-difference 
techniques, each point has its own specific stencil. Although individual stencils may appear to be 
unstable locally, the global method is stable and accurate. Taylor series analysis guarantees that each 
point has a local order property. The use of a specific energy norm to derive the stencils guarantees 
that the resulting global scheme will be stable. 

Another problem with high-order finite-difference techniques, as well as spectral single-domain 
techniques, is their implementation on parallel machines. For conventional finite- difference techniques 
(central or upwind difference), as the order of accuracy increases, the stencil width also increases. 
This results in increased overhead in communicating between processors. Spectral element techniques 
are both efficient and accurate methods for implementation on parallel machines [6). The problem 
is divided into several domains, and each domain is assigned to a processor. Only one point of the 
stencil coexists on multiple processors in spectral element techniques. Information transfer between 
processors is kept to a minimum under these circumstances. 

Cyclo-difference techniques are a combination of high-order finite-difference techniques and spec- 
tral element techniques. They rely on an existing energy-norm proof to establish their stability for 
the hyperbolic system and require no special boundary closure stencils. These techniques can eas- 
ily be split on multi-processor environments. This flexibility results because the discretizations are 
composed of many subelements that each satisfy an order and stability property. The subelements 
are then patched together recursively such that the resulting scheme retains the same stability prop- 
erties. For parallel implementation, they can easily be broken at the patch location, which results in 
a minimum of communication for an arbitrarily high-order scheme. 

SUMMATION-BY-PARTS ENERGY NORM 

Stability of Continuous System 

As shown in reference [5], the summation-by-parts energy norm mimics, at the semidiscrete level, 
the continuous behavior of the principle of the conservation of energy. Because the entire foundation 
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of the cyclo-difference methodology is based on this norm, a complete derivation of the conservation 
of energy principle in the continuous and the discrete case is presented. The model problem is the 
hyperbolic equation defined by 


dU dU_ 
dt dx 


0 


0 <x < 1, t > 0 


( 1 ) 


u(o ,t)=m t> o 


( 2 ) 


U(x, 0) = 4>(x) 0 < x < 1 


( 3 ) 


We begin by defining an energy as E(t ) = U 2 , t > 0. If the energy is differentiated with respect to 
time and the values of U t from equation (1) are substituted, then integration over the domain yields 



t > 0 


( 4 ) 


The definite integral in equation (4) is performed to yield 

E t (i) = -[t/ 2 (M) - U 2 (0,t)] (5) 


and the boundary conditions are substituted from equation (2) to yield 

E t (t) = ~[U 2 (l,t) (6) 

If certain properties on the boundary condition f(t) are assumed, the equation is dissipative (decays 
energy) for all time. For example, if f(t) = 0, then the system energy is uniformly dissipative. 


Stability of Discrete System 

Discrete spatial operators that satisfy very specific properties are shown to be stable in a manner 
analogous to that used in the previous proof of stability. These operators satisfy the summation- 
by-parts energy norm. For example, given the scalar hyperbolic equation U t + U x = 0, a general 
semidiscretization can be written as U t + A* U = 0, where A * is the spatial discretization matrix that 
is presumably consistent to some order. This matrix A* can, in general, be decomposed into the form 

A 9 = P- 1 Q (7) 

This decomposition is in general not unique. If a decomposition can be found such that 


1. Symmetric P : (P = P T ) ( = p jti ) 

2. Positive definite P : ( W T P W > 0) 
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3. Nearly skew-symmetric Q : ( qij = —q h i + 2^,i Sij q lt i + 2 £,,;v £jvj tyv.w) 

4. > 0 and ? u = -qN,N 

then the discretization matrix A* automatically satisfies the summation-by-parts energy norm. (See 
reference [5].) 

To illustrate that this stability property results directly from the form of the matrices P and Q , 
a proof is presented for the semidiscrete form defined by equations (1) and (2). Note that the spatial 
discretization operator can be written in the form 

PU X = QU; U x = P~ l QU (8) 

where P -1 exists and U is the vector of discrete values (f/i, U 2 , U 3 , . . . Un- 2, Un- j, Un) • The semidis- 
crete version of equation (1) becomes 

P^- + QU = 0 t> 0 (9) 

at 

We define the discrete energy as 

E(t) = (U T PU) t > 0 (10) 

where P must be positive definite to ensure that E ( t ) is a strictly positive number. Equation (10) 

t > 0 _ (11) 


0 ( 12 ) 

The semidiscrete expression [equation (9)] is substituted into equation (12) to yield 

Et(t) = 2 [-U t QU] t> 0 (13) 

By using the matrix Q and the relationship between the values q 0 ,o and qx,N, one obtains an energy 
of the form 

£,(() = -2 - Ul) (14) 

The boundary condition defined in equation (2) is substituted to yield 

E t (t) = -2 q N ,N{U 2 N - P(t)) (15) 

Note that the time rate of change of the discrete energy defined in equation (15) is identical in form 
(to within a positive constant) to that of the continuous case equation (5). 
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is differentiated with respect to time to yield the expression 

Because P is symmetric (P = P T ), equation (11) becomes 

E,(t) = ![(/ r ff| (> 


CYCLO-DIFFERENCING 


Conventional central- and upwind-difference techniques use one stencil for the inner portion of 
the spatial domain and auxiliary formulas at the boundaries such that the resulting scheme is stable. 
Spectral element techniques use orthogonal basis functions to define local elements and then connect 
them with various methods that range from spectral patching to flux balances [7]. The simplicity 
and stability of the cyclo-differencing relies on a very specific property of the summation-by-parts 
energy norm. Assume for a particular set of discrete points Xj j = 1, N that a stable discretization 
has been isolated that satisfies all criteria of the summation-by-parts energy norm. The resulting 
semidiscretization written in matrix form is P U t = Q U, where 


Pl.l 

Pi, 2 

. • Pl,N-l 

Pl,N 



P 1,2 

P2,2 

• • P2,N-l 

P2,N 


■ dUi/dt 


. 


dUi/dt 

• 

■ 

• 

■ 

; Ut = 

dUw-\/dt 

Pl,N-l 

P\,N-2 • 

. • PN-l,N-l 

PAT-1, AT 


dU N /dt 

Pl,N 

P2,N ■ 

. • PN-1,N 

PAT, AT j 




-9at,at 9i,2 

—91,2 0 


Q = 


_i_ 

A 


-qi,N-2 

-q\,N - 92 , JV 


qi,N-\ qi,N 

92.AT-1 92, AT 


0 9AT-1.AT 
-qN-i,N qN,N 


' Ui 

u 2 

U N - 1 

u N 


Note that P is symmetric (and positive definite) and Q is skew symmetric except for the elements 
q hl and q N ,N, where q ti i = - qu,N- Also note that the grid-spacing A has been factored out of the 
matrix Q. No attempt has been made to specify the value of the constant N. A cyclo-difference 
scheme is constructed by recursively patching the stable base schemes together, which is illustrated 
by the following example. Assume that the discretization involved 2N - 1 uniformly distributed 
points instead of N points. We can use the properties of the matrices P and Q (of dimension N) 
to define a discretization over the 2N - 1 points that satisfies all criteria of the summation-by-parts 
energy norm, and is therefore stable for any systein of hyperbolic equations. We define this new 
semidiscretization on the equation U t + U x = 0 as PU t = Q U and construct it as 
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where U = [t/i, U 2 , • . • Un-i, Un, Un+ 1 , . . . U 2 N -21 ^ 2 iv-i] • Note that the qw,N element is precisely 
zero because the contributions from qi it and qjv,jv have equal magnitudes but opposite signs. These 
new matrices P and Q satisfy the summation-by-parts energy norm for precisely the same reasons as 
the original matrices P and Q. If the matrices are assembled in this manner, then the matrix P is 
symmetric and the matrix Q is skew symmetric except for the (1,1) and (2 N — 1,2 N — 1) elements; 
these coefficients are again equal in magnitude but opposite in sign. Because the matrix P can be 
decomposed into the summation of two matrices A and A, each of which is positive definite (one, 

A 

but not both, could be positive semidefinite), the resulting matrix P is positive definite. The new 
scheme is therefore stable because the summation-by-parts energy norm is satisfied. In practice, the 
scheme would be implemented as U t = P~ l QU . The inversion of the matrix P could in general be 
quite complicated. 

In short-hand notation, the new matrices P and Q can be written in terms of the original matrices 
P and Q as 
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A 

P = 


P 0 
0 P ; 



A 


Q 0 
0 Q 


This nomenclature is not precisely correct because the new matrix P is (2 N —lx 2 N — 1), and the 
original P is (N x N). More precisely, the last row of the first matrix P and the first row of the 
second matrix P lie on the same row of the new matrix P, with the inevitable intersection of matrices 
at the point pn,n- For the purposes of this work, this nomenclature will not cause any ambiguities. 

Thus far in the derivation, we have assumed that the grid spacing in the first subdomain was 

the same as that in the second domain. In general, this assumption is not necessary, which we will 

demonstrate. Suppose that the first interval is discretized with a grid spacing Aj and the second 

with a spacing A 2 . The resulting semidiscretizations would be P U t = A- Q U and P U t = A Q U, 

respectively. Each respective discretization is multiplied by the appropriate grid spacing to yield 
-£ -# — * — * 

Ai P U t = Q U and A 2 P U t = Q U, respectively. The two subintervals are combined into one to 

yield the matrices P and Q of the form 


P 


A xP 0 
0 A 2 P 



Q 0 
0 Q 


The stability of the resulting scheme is guaranteed by the summation-by-parts energy norm for any 
arbitrary spacing discontinuity. In practice, the scheme would be implemented as U t = P~ l QU. 
The inversion of the matrix P could in general be quite complicated. All information that pertains 
to the discontinuity in spacing is incorporated into the P matrix. 

This procedure of appending new subintervals onto an already existing method can be repeated 
recursively as many times as desired. In fact, even two dissimilar methods (each of which satisfies 
the summation-by-parts energy norm) can be appended to one another on any arbitrary grid-spacing 
interval. One constraint which the resulting cyclo-difference schemes must satisfy, is that the total 
number of grid points must be of the form N t = M (N - 1) +1, where N t is the total number 
of points and M is the number of subintervals that are used. The procedure for grid refinement, 
therefore, involves an increase in the number of subintervals M in the solution. 


HIGH-ORDER CYCLO-DIFFERENCE SCHEMES 


Second Order 

We now present a variety of cyclo- difference schemes of various order and width. In this work, we 
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will concentrate on schemes with uniform grid spacing within each subelement (but not necessarily 
the same from element to element). The first scheme of practical interest is the second-order scheme 
defined on a subelement of three grid points. The scheme represents the optimal second-order scheme 
on three uniformly spaced grid points and, in matrix notation, is given by 


A* 


' -3/2 2 -1/2 ' 

- 1/2 0 1/2 
1/2 -2 3/2 


It is readily shown that A* — P 1 Q, where 



"1/4 

0 

0 ' 


' -3/8 

1/2 

-1/8 ' 

p = 

0 

1 

0 

; Q = 

-1/2 

0 

1/2 


0 

0 

1/4 . 


1/8 

-1/2 

3/8 


This scheme is extended to five grid points to yield P and Q of the form 


'1/4 

0 

0 

0 

0 


' -3/8 

1/2 

-1/8 

0 

0 

0 

1 

0 

0 

0 


-1/2 

0 

1/2 

0 

0 

0 

0 

1/2 

0 

0 

; Q = 

1/8 

-1/2 

0 

1/2 

-1/8 

0 

0 

0 

1 

0 


0 

0 

-1/2 

0 

1/2 

0 

0 

0 

0 

1/4. 


0 

0 

1/8 

-1/2 

3/8 


is 


A* = 


matrix t 

>, P- 

-1 is easily found; the 

' -3/2 

2 

-1/2 

0 

0 

-1/2 

0 

1/2 

0 

0 

1/4 

-1 

0 

1 

-1/4 

0 

0 

-1/2 

0 

1/2 

0 

0 

1/2 

-2 

3/2 


The procedure is extended recursively to an arbitrary number of subelements to yield a cyclo- 
difference scheme of the form 
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' -3/2 

2 

-1/2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 ' 

-1/2 

0 

1/2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1/4 

-1 

0 

1 

-1/4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1/2 

0 

1/2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1/4 

-1 

0 

1 

-1/4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

. 

. 

. 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

. 

. 

. 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

. 

. 

. 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1/4 

-1 

0 

1 

-1/4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1/2 

0 

1/2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1/4 

-1 

0 

1 

-1/4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1/2 

'0 

1/2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1/2 

-2 

3/2 


This scheme is uniformly second-order accurate and satisfies the summation-by-parts energy norm. 
In addition, the scheme does not rely on auxiliary stencils at the boundaries. The operation count is | 
of the count for the conventional second-order scheme and, as will be shown later, behaves noticeably 
different. 

In general, a closed-form expression for A * will not be available because the inverse of P will not 
be known analytically. Therefore, a banded solver of width 2N — 1 (the number of points in each 
subelement) can be used on the matrix P to efficiently invert the matrix. Although banded solvers 
are efficient in comparison with full solvers, they can not compete with explicit schemes in which no 
numerical inversion of the matrix P must be performed. The previous example demonstrates that 
methods with a diagonal matrix P can be immediately inverted. By only concerning ourselves with 
those numerical methods that possess this property, we are being overly restrictive. In general, if P 
has a first and a last row that consists entirely of zeroes except for the diagonal element, then each 
subelement of the resulting matrix P decouples and the inverse can be performed analytically. The 
resulting scheme has an explicit, not implicit, operation count and can compete with conventional 
finite-difference schemes of comparable spatial accuracy. 

Third Order 

A uniformly third-order scheme can be generated with a minimum of four discrete points. The 
discretization matrix A *, which is third order and occupies four points, can be represented as 

-11/6 3 -3/2 1/3 ' 

-1/3 -1/2 1 -1/6 

1/6 -1 1/2 1/3 

-1/3 3/2 -3 11/6 . 

If A* is decomposed into P~ l Q, the resulting relationships for P and Q are 
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p 


r3 

360 r3— 11 r2-j-21 r 1 
312 

-(54r3-r2+9rl) 

78 

120r3-f5r2+33rl 

312 


360 r3— 11 r2-f 21 rl 
312 

207r3+7r2+15rl 

39 

-(72 r3+55 r2+51 rl) 
312 

— (54 r3— r2-j-9 rl) 
78 


-(54 r3-r2+9 rl) 
78 

-(72 r3+55 r2+51 rl) 
312 

207 r3+7r2+15rl 
39 

360 r3-llr2+21rl 
312 


120 r3+5r24-33rl 
312 

-(54 r3-r2+9rl) 
78 

360 r3— 11 r2+21 rl 
312 


r3 . 


Q = 


-(288r3-r2+9rl) 

-(384 r3+ir2+25 rl) 
104 

24 r3+r2+4rl 

-f576 r3+3 1 7r2+135rl) 
936 


384 r3+3 r2+25rl 
104 

0 


~(576r3+llr2+57rl) 

104 

24 r3-i-r2+4 rl 
13 


-(24 r3+r2+4 rl) 

13 

576 r3+ll r24-57 rl 
104 

0 

-(384 r3+3 r2+25rl) 
104 


576r3+37r2+135rl * 
, 936 

-(24 r3+r2+4rl) 

13 

384r3+3r2+25rl 

104 

288 r3-r2+9 rl 

117 J 


The final criteria to be met is that the matrix P must be positive definite. Even if this scheme 
was stable for arbitrary rl, r2, and r3, the amount of work necessary to invert the matrix P would 
make the scheme prohibitively expensive. Therefore, the free parameters are used to decrease the 
bandwidth of the matrix P. If we set r3 = 1, r2 = 9r3, and rl = — 5r3, then P and Q result of 
the form 


' 1 

1/2 

0 

0 ' 


■ -2 

11/4 

-1 

1/4 ] 

1/2 

5 

-1 

0 

. f) _ 

(— H)/4 

0 

15/4 

-1 

0 

-1 

5 

1/2 


1 

(— 15)/4 

0 

11/4 

0 

0 

1/2 

1 


(-l)/4 

1 

( 11)/ 4 

2 


With Gershgorin’s Theorem, the matrix P is shown to be positive definite. The original matrix 
satisfies the summation-by-parts energy norm and is stable for the hyperbolic system. Unfortunately, 
not enough free parameters exist in the decomposition to make P diagonal. The work involved in using 
a cyclo-difference scheme generated from P and Q would be seven multiplications and additions per 
node: four from the Q matrix and three from the inversion of the tridiagonal matrix P. An operation 
count this high would not make the resulting scheme competitive with other third or higher order 
schemes. 

To obtain high-order schemes in which P can be diagonalized, or at least decoupled from the 
other subelements, nonoptimal schemes must be used. Consider, for example, the family of uniformly 
third-order schemes that can be defined on five uniform points. Because each point allows a new 
degree of freedom, a wide variety of different schemes can be developed. Assuming that the following 
constraints are imposed: (1) uniform third-order accuracy exists from five points,. (2) a P and Q 
exist that satisfy the summation-by-parts energy norm, and (3) a matrix P exists that has a first 
and last row composed entirely of zeros except the diagonal element. Matrices P and Q result of the 
form 
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• 1 

0 

0 

0 

0 ' 



0 

—(13 rl— 352) 

7 rl— 160 

—(7 rl — 160) 

0 



12 

, 6 

12 


p = 

0 

7 rl — 160 
6 

—(25 rl— 592) 
12 

7 rl— 160 
6 

0 

> 


0 

—(7 rl— 160) 

7 rl — 160 

—(13 rl— 352) 

0 



12 

6 

12 



. 0 

0 

0 

0 

1 . 




[ 3 rl— 120 

-(9 rl— 256) 

9 r 1—208 

— (3rl— 64) 

9 rl — 184 


32 

24 

16 

8 

96 


9 rl— 256 
, 24 

0 

24 — rl 

3 rl —64 
3 

— (3rl— 64) 
8 

= 

—(9 rl— 208) 
16 

rl -24 

0 

24 -rl 

9 r 1 — 208 
16 


3 rl —64 
. 8 

—(3 rl— 64) 
3 

rl -24 

0 

—(9 rl — 256) 
24 


— (9 rl — 184) 

3 rl— 64 

—(9 rl— 208) 

9 rl— 256 

-(3 rl — 120) 


L 96 

8 

16 

24 

32 


The characteristic polynomial of the matrix P is 

(A — l) 2 (2 A + rl -32) [12 A 2 + (45 rl - 1104) A + 9rl 2 - 560 rl + 8192] = 0 

for which the roots are A = 1, »=!l and ± g=gSH ±13gg. For values of 

rl < 280 ~ 8v ^ , all eigenvalues are positive, and the method is stable. The resulting matrix A* can 
be written as 



[ 3 r 1 — 1 20 

, 32 > 

— (rl— 48) - 

(9 rl— 256) 
24 

(7 rl — 160) 

9 r 1—208 
16 

2 rl— 48 

—(3 rl— 64) 
8 

— (5 rl— 96) 

9 rl — 184 
96 

rl — 16 

A * = 

6 rl — 192 
1 

, 12 v 
— (rl — 16) 

6 rl— 192 
-2 
3 

5 rl -96 

rl— 32 
0 

— {2 rl— 48) 

6 rl — 192 

2 

3 

7 rl — 160 

6 rl — 192 
-1 

12 

rl— 48 


6 rl — 192 
— (9 rl — 184) 
L 96 

6 rl — 192 
3 rl— 64 
8 

rl-32 
—(9 rl— 208) 
16 

6 rl— 192 
9 rl — 256 
24 

6 rl — 192 
—(3 rl — 120) 

32 

Note that for a value of rl = ^ the resulting scheme is 




P = 

' 1 0 
0 32/7 
0 0 
0 0 
0 0 

0 0 0 ' 

0 0 0 

12/7 0 0 ; 

0 32/7 0 

0 0 1 


Q = 

' (— 45)/28 
(— 44)/21 
1/7 

4/7 i 

. ( — 19)/84 

44/21 

0 

(— 8)/7 
(— 32)/21 
4/7 

(~l)/7 

8/7 

0 

( 8)/7 
1/7 ( 

(— 4)/7 
32/21 
8/7 
0 

— 44)/21 

19/84 ' 
(— 4)/7 
(-l)/7 ; 
44/21 
45/28 
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' (— 45)/28 

44/21 

(~ l )/7 

(" 4)/7 

19/84 ‘ 

( — ll)/24 

0 

1/4 

1/3 

( l )/8 

1/12 

(— 2)/3 

0 

2/3 

(-m 2 

1/8 

(-l)/3 

(- l )/4 

0 

11/24 

(— 19)/84 

4/7 

1/7 

( 44)/21 

45/28 


Other third-order subelements exist for a uniform grid, but require a larger stencil. Their operation 
count is necessarily larger than those already presented and will not be pursued in this work. 


Fourth Order 

The optimal fourth-order schemes defined on five grid points produce a subelement A * that can be 
decomposed into P and Q of the form 


1 1 

2/3 

0 0 

0 


2/3 

20/3 ( 

-10)/3 4/3 

0 


0 (- 

-10)/3 

32/3 (— 10)/3 0 


0 

4/3 ( 

— 10)/3 20/3 

2/3 


0 

0 

0 2/3 

1 


r (— 9)/4 

31/9 • 

-2 1 ( 

-7)/36 ' 


( 31 )/9 0 

6 (— 32)/9 

1 


2 

-6 

0 6 

-2 

1 

-1 

32/9 • 

-6 0 

31/9 


7/36 

-1 

2 ( — 31)/9 

9/4 


— 25)/12 

4 

— 3 4/3 

(— 1)/4 

( 1 )/4 

(-5)/6 

3/2 (— 1)/2 1/12 

1/12 

(-2)/3 

0 2/3 

(— 1)/12 

( 1)/12 

1/2 

( — 3)/2 5/6 

1/4 

1/4 

(-4)/3 

3 -4 

25/12 


Gershgorin’s theorem can be used to prove that the matrix P is positive definite. The scheme satisfies 
all criteria of the summation-by-parts energy norm and is, therefore, stable for hyperbolic systems. 
Note that the resulting scheme is pentadiagonal in the matrix P and would require a banded solver in 
the cyclo-difference mode to invert. The operation count of just the P matrix inversion would be 5 N, 
which is not competitive with other explicit formulations. As the order of accuracy increases with 
optimal formulations, the bandwidth of the matrix P will also increase, which makes these schemes 
impractical. In addition, for N sufficiently large, the resulting schemes are unstable and cannot be 
used. 
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A subelement of uniformly fourth-order schemes that is stable and explicit can be derived by 
considering six uniformly distributed discrete points. One additional degree of freedom from each 
point enables all of the constraints to be met. As with the third-order subelements, a matrix P that 
decouples is sought. Matrices P and Q that satisfy these constraints are 


r i 


0 


0 


0 


0 


0 n 


0 

0 

0 

0 


—(2797 rl -40008700) 

37TOT5 

287 rl -2348300 

4 163636 ' 

-(249 rl — 1862000) 

* TSSU7* 

287 rl -2348300 

— mm — 


287 rl -2348300 
183028 

-(2999 rl -20919400) 

irmb 

557 rl —3154500 

163626 

-(249 rl -1862000) 

— mm 


-(249 rl -1862000) 
— t 

557 rl -3154500 
OT5 

-(2999 rl -20919400) 

913146 

287 rl -2348300 

— mm — 


0 


0 


0 


0 


287 rl -2348300 

— mm — 0 

-(249 rl -1882000) 

— mm 0 

287r fe«r 30Q o 

—(2797 rl -40008700) 

2745426 0 

0 1 


r 72 rl -2655360 

1)43955 

864 rl —14247875 

-(864 rl -14247875) 
2745426 

864 rl -7384325 

T t mm — 

-(1464 rl — 1 1209525) 

-(288 rl — 1698825) 
457570 

2028 rl -6385925 

432 rl —1976275 

1375716 

-(952 rl —2851075) 

-(864 rl -3266195) 
‘13727166 

432 rl —1976275 

274542b 


1372716 

1372710 

9'irHO 

1372710 

-(864 rl -7384325) 

1464 rl -11209525 


-(882 rl -2557325) 

2028 rl -6385925 

-(288 rl —1698825) 

1372716 

1372716 


666355 

1372716 

457570 

288 rl -1898825 

-(2028 rl -6385925) 

882 rl -2557325 

n 

-(1464 rl — 11209525) 

864 rl -7384325 

457570 

1372716 

886355 


“ *“ 1373710 

1372716 

-(432 rl— 1976275) 

952 rl —2851075 

-(2028 rl -6385925) 

1464 rl -11209525 


-(864 rl —14247875) 

1372716 

~ 915140 

1372710 

“ i372YiiT 

0 

J745"426 

864 rl -3266195 

-(432 rl —1976275) 

288 rl— 1698825 

-(864 t 1—7384325) 

864 rl —14247875 

-(72 rl -2655360) 

L 13727100 

1 372710 

457570 

1372710 ~~ 

2745420 ““ 

1143928 


The roots of the characteristic polynomial of the matrix P are A = 1, 167 ^ 5 1 ^q 42325 ± 

5 \/5 n/X0T7i 3 -2874060 rl +20964334625 J -5396 rl +40456475 i 5 ^1076752 rl?- 13981396400 r 1+483895685 78125 p 

457570 ’ 13727 10 - 1 - 137271Q * rOP 

values of rl < 36 ^[g| 75 — l 25 ( a numerical value of about 5733), all eigenvalues are 
positive, and the method is stable. 


Fifth Order 

Optimal subelements of fifth-order accuracy are not pursued in this work. Instead, a uniformly fifth- 
order accurate scheme defined on seven evenly spaced points is derived. The matrices P and Q are 
defined by 


1 0 

n 1 06 9 r 1 4292170384 
25286256 


-(89 rl + 7400232) 

0 - 1 — 1652564 

A 863 rl + 44908416 
0 8478755 

„ -(863 rl + 44908416) 

0 138*5178 

o 2^^232 

0 0 


-(89 rl + 7400232) 

— 1653594 L 

2627 rl + 11311560 

12653128 

-(1111 rl —41866272) 

4214376 

809 rl -13279392 

4214376 

-(863 rl + 44908416) 

12643)38 


863 rl+44908416 
— — 

-(1111 rl —41866272) 

4214376 


1489 r 1-72846720 
4214376 

-(1111 rl —41866272) 

4214376 


863 rl+44908416 

— mV — “ 


28752 


-(863 rl+44908416) 

12643128 — 

809 rl — 13279392 

4214376 

-(1111 rl —41866272) 

4214376 * 

2627 rl + 11311560 
12643128 

-(89 rl + 7400232) 

* 1653594 L 


89 rl + 7400232 

4-2)4376 


-(863 rl+4490B416) 
12643128 


863 rl+44908416 
PT25 


8428752 

-(89 rl + 7400232) 

msm — 

1069 rlj-2921 70384 


0 
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Q - 


-(325 rl+43031322) 
' 10535340 L 


125 rl+II423502 
4214376 

-(125 rl+887562) 
— 3T507BS 


125 rl —4380408 
471*375 

-(25 rl -1508238} 
— 2107188 

. 125 rl —9648378 
63215640 


125 rl + 43031322 
10535940 


-(1255 rl -486938) 
25286256 

185 rl -28229184 
2107188 

-(715 rl —127212408) 
8*28752 


1475 rl -205934976 

3*60 7 820 

-(25 rl -1508238) 
— 2'107I88 


-(125 rl + 1 1423502) 
4214376 u 


1255 rl -486936 
25286256 


0 

-(305 rl -88569864) 
4214376 L 

680 rl -204436503 
8321564 

-(715 rl -127212408) 
8*28752 

125 rl -4380408 
4214376 



125 rl + 887562 
" 3160782 

-(125 rl —4380408) 
4214376 L 

25 rl — 1508238 
2107188 

-(125 rl -9648378) "I 
63215640 


-(185 rl -28229184) 

715 rl — 127212408 

-(1475 rl -205934976) 

25 rl -1508238 


2107188 

8428752 

31807820 

2107188 


305 rl —88569864 

4214378 

-(680 rl -204436503) 
833TO4 L 

715 rl —127212408 
8478752 

-(125 rl -4380408) 
*27*576 L 


... o 

305 rl -88569864 

-(185 rl— 28229184) 

125 rl+887562 



4214376 “ 

2107188 

3160782 


-(305 rl -88569864) 

o 

1255 rl —486936 

-(125 rl + 1 1423502) 


— 4214376 ‘ 


25288256 

4214376 


185 rl —28229184 

-(1255 rl -486936) 

A 

125 rl + 43031322 


2107188 “ 

25286256 

u 

10533940 


-(! 25 rl+887562) 

125 rl + 11423502 

-(125 rl+43031322) 

25 rl + 26938800 

- 

3160782 

*21*376 

““ 10535940 

' — T23f3T28 J 


The matrix P is positive definite for values of rl > 172357 (determined numerically). This 
uniformly fifth-order scheme satisfies the summation-by-parts energy norm and can be used as the 
subelement for generating a cyclo-difference scheme. The resulting scheme will be explicit in nature 
because the matrix P can be inverted analytically. 

This basic procedure can be used to generate schemes of arbitrarily high order, although none 
with an accuracy of greater the fifth order was generated in this work. 


STABILITY OF THE CYCLO-DIFFERENCE SCHEMES 


The summation-by-parts energy norm was used to develop the cyclo-difference schemes in a 
semidiscrete context. The theoretical CFL for various Runge-Kutta (R-K) schemes must still be 
determined. Because each point in a cyclo-difference scheme uses a different stencil, the use of 
Fourier techniques to obtain a CFL is not applicable. A numerically determined eigenvalue spectrum 
provides a practical means of obtaining the CFL of the various schemes. 

A necessary condition for Lax stability of a semidiscretization is that the eigenvalues of the spatial 
discretization operator (scaled by At), lie within O(Atf) of the stability region of the time integration 
formula as At — ► 0. This condition allows for exponential growth of the solution in time, but at a 
rate that is independent of the grid density, which guarantees convergence as Ax — » 0. For time 
asymptotic stability, all eigenvalues of the spatial discretization operator (scaled by At) should lie 
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within the stability region of the time integration formula for all At. This condition guarantees 
that the solution in time will remain bounded if the solution to the original governing equation is 
bounded. The determination the CFL of a discretization is thus reduced to solving for At such that 
the resulting stability region of the time integration formula encompasses all spatial eigenvalues. 

Figure 1 shows the eigenvalue spectrum of the cyclo-difference schemes as determined from a 
numerical eigenvalue determination. For each case, the number of grid points is 61; this number 
satisfies the grid constraints for all cyclo-difference schemes presented thus far. Note that the structure 
of the spectrum is not continuous, but seems to cluster into specific portions of the complex plane as 
for conventional finite-difference techniques. (See reference [2].) In addition, this clustering becomes 
more pronounced as the order of accuracy increases. 

Table I shows the CFLs of the various schemes determined from a numerical eigenvalue deter- 
mination. In all cases, the grid used contained 61 points. Very little sensitivity to grid density was 
observed between 31 and 61 grid points. The four algorithms used in the study were : 1) cyc23, 2) 
cyc35, 3) cyc46, and 4) cyc57. The first number indicates the formal accuracy of the scheme; the 
last number is the number of grid points occupied by the subelement stencil. The cyclo- difference 
schemes are then generated by recursively appending the subelements. Note that each scheme can 
only run on a grid of M(N — 1) -f 1 points, where N is the element size and M is the number of 
elements. 


(R-K) 

CFL 

CFL 

CFL 

CFL 

Order 

Cyc23 

Cyc35 

Cyc46 

Cyc57 

2nd 

3rd 

4th 

Unstable 

1.16 

1.89 

Unstable 

1.23 

2.01 

Unstable 

1.08 

1.77 

Unstable 

0.87 

1.42 


Table I. CFL of cyclo- difference schemes determined from eigenvalue determination. 


Note that in each spatial operator, the fourth-order R-K is more efficient (in terms of CPU time, 
rather than storage) than the third-order R-K to advance the solution in time. 

Interestingly, equation (14) was obtained without imposing the boundary conditions. Reference [5] 
demonstrated that the method of imposing physical boundary conditions is important to the overall 
stability of the method. Two counter examples showed that the property of time asymptotic stability 
could not in general be guaranteed for discretizations that satisfy the summation-by-parts norm. 

The underlying reason for growth in time is the effect of the boundary operator on the structure 
of the norm matrix P -1 . Recall the semidiscrete form of equation (1): Ut = P _1 Q- If the boundary 

15 



. , , ■ f,, rTn tt _ n p- 1 n results, for which DP 1 may not be a 

operator D is imposed, the semidiscrete form U t - u r v resui , * 

norm. In the scalar case, this result is not a concern if the matrix P is a restarted full norm (Strand 

I IP or if D P - 1 = P" 1 D. A restricted full norm is defined when the diagonal is the only nonzero 

element in the first (or last) row and column of the matrix P. The cyclo-difference schemes satisfy 

precisely the definition of the restricted full norm and, therefore, are automatically stable for the 

scalar wave equation. 

Unfortunately, even for cases where P is a restricted full norm, stability cannot be generalised 
to the case of a general hyperbolic system. The SAT method proposed in reference [5] can be used 
to guarantee stability for the hyperbolic system. The SAT method involves the indirect imposition 
of the boundary conditions, which is accomplished by adding a term to the derivative operator that 
is proportional to the difference between the discrete value U„ and the boundary term /((). Rather 
than directly satisfying the boundary condition by imposing U N = f(t), the SAT method solves a 
derivative equation throughout the entire domain including the boundary points. 

ORDER PROPERTIES OF CYCLO-DIFFERENCE SCHEMES 


Time Dependent 

Several test problems (both steady and unsteady) are used to establish the accuracy of the cyclo- 
difference schemes. Because Taylor series analysis was used in all cases to derive the schemes o a 
particular order, these schemes are expected to behave to at least the order of the local truncation 
error. Consider the method-of-lines approximation to the scalar wave equation 


d JL + d JL = o 

dt dx 


-1 <x < 1, <>0 


1 < x < 1, * >0 


(16) 

(17) 


U(t, —1) = sin 2 tt( — 1 — <); t/(0,x) = sin2?rx 
with the exact solutions given by 

U(t,x) = sin2ir(z — f), -1 < * < 1, > > « < 18 ) 

The spatial discretization is accomplished by the new cyclodifference schemes of various order; time 
was advanced with a four-stage R-K time-advancement scheme (formal nonlinear accuracy of fourth 
order) with a CFL of 0.25 to a time level of 25. Further temporal refinement showed no improvement 

in the solution accuracy. 

Uniform Mesh . . , „ 

We begin with a discretization on a uniform mesh. Table II shows the results from a grid-refinement 

study performed with each algorithm. 
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Grid 


log 




Cyc23 

Cyc35 

Cyc46 

Cyc57 

16 



-1.347 


17 


-0.7378 



19 




-1.530 

21 

-1.662 

-0.6959 

-1.740 


25 


-1.625 


-2.059 

29 


-2.066 



31 

-2.083 


-2.276 

-2.760 

33 


-2.218 



37 




-3.273 

41 


-2.587 

-2.722 


61 

-2.710 

-3.521 

-3.736 

-4.620 

81 


-3.914 

-4.249 


91 



-4.455 

-5.678 

101 


-4.248 

-4.638 


121 

-3.316 

-4.547 

-4.956 

-6.435 


Table II. Accuracy of new cyclo-difference schemes as function of grid density. 


Once a scheme has achieved a certain minimum grid density it will exhibit an order property; by 
doubling the mesh, the error will decrease by a factor of 2*, where k is the order of the scheme. The 
slopes for each scheme, as determined between grid densities of 31 (37 for cyc57) and 121 points, 
were -2.05, -4.05, -4.45, and -6.10, respectively. 

Note that the apparent accuracy of some of the cyclo-difference schemes is higher than the pre- 
dicted local truncation error. This result is more apparent in the schemes with an odd order, namely 
the cyc35 and the cyc57. The cyc46 does not obtain an accuracy that is significantly greater than 
the theoretical accuracy. Borrowing from the finite-element terminology, this increased convergence 
rate shall be referred to as “superconvergence.” 

Discontinuous Mesh 

We noted earlier in the derivation of the underlying principles of cyclo-differencing that two subinter- 
vals of unequal spacing could be joined at an interface without destroying the accuracy or stability 
of the formulation. The discretization matrices P and Q were defined by 


P = 


A X P 0 
0 A 2 P 



Q 0 
0 Q 


Table 111(a) shows the results of a grid-refinement study for which the ratio of grid spacings between 
subintervals was not one ^ 1). The model problem was that used in the previous test case 
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(equations (16) to (18)), and the spatial discretization was the cyc35 algorithm. Time advancement 
was as previously described. In each case, the spatial domain was divided into two regions, and one 
half of the total number of points was distributed uniformly throughout each domain. This gave rise 
to grid spacings A1 and A2 in each domain, respectively. Table 111(a) shows the logarithm of the 
error for each discretization as a function of grid density, for a variety of spacing ratios p = ^ • In 
all cases, the finest concentration of mesh points occurred at the outflow boundary. 


Grid 


logI 2 


HNH 


EHI 

p = 3/2 

EBI 

EBQ 

41 

-2.587 

-2.469 

-1.864 

-1.849 

81 

-3.914 

-3.715 

-3.197 

-3.123 

121 

-4.547 

-4.521 

-3.974 

-3.915 

161 

-5.046 

-5.052 

-4.559 

-4.315 


Table 111(a). Accuracy of cyc35-difference scheme for discontinuous mesh spacing as function of 

grid density. 


The magnitude of the error increases as the mesh discontinuity increases. This is because of the 
increased effective grid density that results from clustering the mesh points near the outflow boundary. 
The scheme still behaves with a fourth-order accuracy on this problem. Note that the slope of error 
decay for the p = 5 case is -3.96 between N - 81 and N = 161 and that the cyc35 scheme is still 
superconvergent . 


Table 111(b) shows a similar comparison with the other cyclo-difference algorithms; all cases were 
run with a grid ratio of p = 5. 


Grid 


logZ -2 



cyc35 

cyc46 

cyc57 

41 

61 

-1.849 

-2.317 

-3.251 

81 

97 

-3.123 

-3.279 

-4.449 

121 

-3.915 

-3.983 

-5.004 

161 

-4.315 

-4.476 



Table 111(b). Accuracy of cyclo-difference schemes for discontinuous mesh spacing (p = 5) as 

function of grid density. 


The slope of the error decay in the cyc57 scheme between points 61 and 121 is -5.82. The odd-order 
schemes, even in this discontinuous case, appear to be superconvergent. 
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Steady State 


The second test problem is the solution of the flow through a supersonic nozzle. The governing 
equations are the quasi-one-dimensional Euler equations. For this problem, an exact steady-state 
solution exists and can be used to compare the accuracy of the new cyclo-difference methods. The 
governing equations are 


JL(M) + |;(,uA) = 0 

+ p) A] = p |i 

^{pe t A) + -^[(pe t +p)uA] = 0 (19) 

where A = A(x) and e t = Cv T + y . Boundary conditions are imposed on the inflow plane for all 
three equations, and A(x) is prescribed such that the flow remains supersonic throughout the entire 
domain. A four- and a five-stage R-K were used to time advance the solution to the steady state 
(machine precision of 10 -13 ). Residual smoothing was used to accelerate the convergence rate for the 
various schemes. Table IV shows a comparison of the Li error that resulted from each of the four 
cyclo-difference schemes on various grids. 


Grid 


log £ 2 




Cyc23 

Cyc35 

Cyc46 

Cyc57 

121 

-4.122 

-5.450 

-6.410 

-6.851 

181 

-4.474 

-6.153 

-7.178 

-7.916 

241 

-4.723 

-6.653 

-7.703 

-8.666 


Table IV. Accuracy of Cyclo-difference schemes as function of grid density for one-dimensional 

nozzle flow. 


Note that by doubling the mesh, an error decay is produced with a slope of -2.00, -4.00, -4.29, and 
-6.00, respectively. These slopes agree with those obtained in the time-dependent case for simple 
linear advection. Again we see that the odd-ordered schemes are superconvergent. 

Note that conventional second-order methods (as well as higher order central methods) with 
suitable boundary conditions will not converge to steady state for this and many other practical flow 
problems. The residual decreases only one order and then remains at this point indefinitely. This 
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non convergence is because the interior scheme is entirely dispersive; thus, the only dissipation in the 
spatial scheme comes from the boundary closure terms, which generally are not sufficient to damp 
the odd-even ( tt 7r” in Fourier space ) mode that can develop under non-linear circumstances. Higher 
order damping is explicitly added to these central-difference schemes to ensure that a steady-state 
solution can be found. The cyclo- difference schemes, by virtue of their cellular construction, are 
stable for the ir mode and can be used without additional damping. In practice, additional damping 
makes the scheme converge more rapidly as expected. 

CONCLUSIONS 

A new methodology referred to as “cyclo-differencing” is presented for defining stable high-order 
finite-difference schemes. They are similar to spectral element techniques of some arbitrary order, 
and are ideally suited for implementation on parallel machines. Unlike spectral element techniques, 
their existence does not rely on orthogonal polynomials or nonuniform grids. In principle, they can 
be devised of any order accuracy, although only schemes of sixth-order accuracy are pursued in this 
work. These new techniques rely on the summation-by-parts energy norm to establish formal stability 
for the scalar hyperbolic case. In addition, if the newly devised SAT method for imposing boundary 
conditions is used in conjunction with cyclo-differencing, then the resulting numerical method is 
formally stable for the hyperbolic system. 

The cyclo- difference techniques are similar to central-difference techniques in that they are stable 
for right- or left-running waves, with appropriate placement of the physical boundary conditions. A 
decided advantage is that no numerical boundary conditions are required near the walls. Thus, high 
order accuracy is assured throughout the entire domain. In addition, the cyclo- difference schemes can 
be patched together across arbitrary grid discontinuities, and still retain their accuracy and stability. 

A series of test problems are used to demonstrate the efficacy of the cyclo-difference methodology. 
The scalar advection equation is used to show the formal stability and accuracy of the second- 
to fifth-order cyclo-difference schemes. For the odd-order cyclo-difference schemes, the property of 
superconvergence is observed; specifically, these schemes converge at a rate one order higher than their 
theoretical accuracy on both uniform and discontinuous grids. A one-dimensional nozzle problem 
is used to demonstrate the robustness of the cyclo-difference techniques. Steady-state solutions, 
consistent with the order property of the spatial operator, are obtained without the addition of 
artificial damping to the formulations. Finally, the viscous Burgers equation is solved to demonstrate 
the use of the cyclo-difference technique for a nonlinear parabolic problem. Again, robust and accurate 
solutions are obtained in all cases. 
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APPENDIX A 


Viscous Stability 

Assume P = QU, where P and Q are matrix operators and U is a vector of discrete values. The 
proof of stability for the cyclo- difference scheme developed in this work relies on very specific forms 
for the matrices P and Q. Specifically, P = P T and is positive definite; Q = Q, m + Q,k, where the 
only nonzero elements of Q a k are 90,0 and qN,N, and <7o,o = —Qn.n = —a; (a > 0). This form of the 
derivative operator leads directly to a stable second derivative operator. We begin with a derivation 
of the continuous energy for the one-dimensional heat equation 


dU d 2 U 
dt ~ Qd dx 2 


0 <x < 1, t >0 


U(0,t) = f(ty, U(l,t)=g(t) 


U(x,0) = ip(x) 0 < a: < 1 (A. 1) 

Note that boundary conditions based on the derivatives at x = 0 or x = 1 could have been imposed. 
An energy (defined as \ U 2 ) is formed by multiplying equation (A. 1) by U. Integration over the 
domain results in 

E t W = Jj~ ad U dX t ~° (A ‘ 2) 

Integration by parts yields an expression of the form 

E t (t) = a d [U(l,t)U x (l,t) - U(0,t)U x (0,t)] - J a d (^) 2 dx t> 0 (A. 3) 

The energy takes the form of a negative definite quantity plus the boundary data that involve the 
function and its derivatives at the boundaries. 

If the second derivative in equation (A. 1) is formed by twice operating with the first derivative 
operator, the semidiscrete form of equation (A. 1) becomes 

^ = a d P~'QP- x QU t > 0 (A. 4) 

If P is symmetric ( P = P T ) and is positive definite, then P~ x is symmetric and positive definite. 
Similarly, because Q = Q aym + Q a k-, we have Q = 2 Q sym - Q T . By operating on equation (A. 4) 
from the left by U T P and using the relationships between Q and Q T , we obtain 

E t {t) = a d \2U T Q T tym P~ l QU - U T Q T P~ X QU] t> 0 (A. 5) 
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where E t is the time rate of change of the discrete energy defined by E(t) = U T P U. In defining 
V = QU, the second term on the right side of equation (A. 5) is negative definite. Because of the 
sparseness of the matrix Q^ ym , the first term reduces to a [U(l) fj(l) — U(0) §^((J) ] and equation 
(A. 5) becomes 


E t (t) = aj 


“( 




U(0) 


» - 


V T p- 1 V 


t> 0 


(A. 6) 


To within a positive constant, this expression for the discrete energy is identical to that obtained from 
the continuous case. Therefore, the discrete energy will behave like the continuous energy. However, 
note that this analysis is linear and does not guarantee stability in some problems of practical interest. 
For example, the same argument could have been used to demonstrate the stability of forming U 2x 
from two central-difference first derivatives. In fact, if U 2x is formed in this manner, no damping 
would result for the 7 r mode in Fourier space, which ultimately would result in the growth of the 
odd-even mode in physical space. 


The cyclo-difference schemes do not appear to be as susceptible to the odd-even mode instability 
as the conventional central-difference schemes because a different stencil is used at each point. For 
example, consider the cyc35 scheme with the parameter rl = 16. The resulting subelement A * = 
P- 1 Q is 


' -9/4 

14/3 

-4 

2 

-5/12 

-1/3 

-1/2 

1 

-1/6 

0 

1/12 

-2/3 

0 

2/3 

-1/12 

0 

1/6 

-1 

1/2 

1/3 

5/12 

-2 

4 

-14/3 

9/4 


By forming the viscous derivative by the sequential operation of the first derivative operator A v = 
A * A % we obtain 



3 -9 

10 

-5 

1 

1 -2 

1 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

1 

-2 

l 

1 -5 

10 

-9 

3 


The interior of this matrix is identical to the conventional second-order viscous derivative and is 
not susceptible to the odd-even mode. Truncation analysis shows the resulting viscous matrix to be 
uniformly second order. 
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To test the convergence rate of the viscous terms and the overall influence of the odd-even mode, 
the viscous Burgers equation is used. The equation is defined by U t + (0.5 U 2 ) x = fi U xx , with 
boundary conditions 17(0,/) = U 0 and 17(1,/) = 0. The exact steady-state solution is given by 


U(x) = U 0 U 


1 — exp U Re (x — 1) 
1 + exp 17 Re (x — 1) 


where Re = Uo/fi and U is the solution of the equation 

^ = exp (— U Re) 

U + 1 ' 

and can be used to determine the error on a particular grid for this problem. 


Table A. I shows the results of a grid- refinement study on the viscous Burgers equation. 


Grid 

fourth 

Cyc35 

Cyc46 

Cyc57 

51 

53 

55 

-2.009 

-3.405 

-3.547 

-2.975 

101 

103 

-3.408 

-4.663 

-4.496 

-4.719 

201 

205 

-4.879 

-5.882 

-6.233 

-6.807 

401 

403 

-6.369 

-7.087 

-7.561 

-8.596 


Table A. I: Error from various cyclo- difference schemes on the steady-state Burgers equation. 

The log 10 of the error are plotted as a function of the grid density. The convergence rate for 
the fourth, Cyc35, Cyc46, and Cyc57 methods are 4.9, 4.0, 5.1, and 6.5, respectively, as determined 
between the 101 and 401 points. The viscosity was n = 0.04 which results in Re = 25. 

APPENDIX B 


Periodic Stability 

At least two cycles of the fundamental subelement are required to define a cyclo-difference scheme 
(which includes two boundaries and one patch). Schemes of greater grid density can be constructed 
by cyclically patching an arbitrary number of subelements together. The cyclo- difference schemes 
can be used, however, on a periodic domain. By construction, we will show how to generate a stable 
periodic scheme from any of the cyclo- difference schemes. 
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The periodic assumption is implemented by first requiring that grid points 1 and N are equivalent. 
Hence, the last row and column from the cyclo-difference scheme can be eliminated. (With this 
requirement, the minimum number of subelements required for the periodic case is three.) If points 
1 and N are equivalent, then the stencils that require grid point N are “wrapped around” to point 1. 
Similarly, the stencil at point 1 is replaced with the interface stencil symmetrically relates points on 
either side of grid point 1. The resulting cyclo-difference stencil is now periodic, and each subelement 
is indistinguishable in terms of position. Unlike conventional central-difference schemes, the resulting 
stencil is not skew symmetric. The eigenvalues of the stencil are all on the imaginary axis because 
A p = P~ l Q P , and Q v is entirely skew symmetric. 

The periodic versions of the cyc23 and cyc35 schemes are presented to illustrate this procedure. 
The periodic cyc23 scheme with three subelements (seven points) and truncated to six points, can 
be written as 


A p 


' 0 

1 

-1/4 

0 

1/4 

-1 ' 

-1/2 

0 

1/2 

0 

0 

0 

1/4 

-1 

0 

1 

-1/4 

0 

0 

0 

-1/2 

0 

1/2 

0 

-1/4 

0 

1/4 

-1 

0 

1 

1/2 

0 

0 

0 

-1/2 

0 


The matrix A p = P~ x Q p , where P is the diagonal matrix characterized by [1/2, 1, 1/2, 1, 1/2, 1] on 
the main diagonal; Q p is entirely skew symmetric. By definition, we know that a skew symmetric 
matrix has eigenvalues on the imaginary axis. The semidiscrete energy of the system defined by 
U T P p U will be unchanged for all time with this discretization. 


A comparison of the eigenvalues for this discretization with those of the conventional second- 
order periodic central-difference stencil on six points is interesting. The characteristic polynomial 
for the matrix A p is A 2 (16 A 4 + 51 A 2 + 36) = 0, which results in the eigenvalues A = 0 and 
A = ± ^3/32 \j\l ± \/33i for which the maximum eigenvalue is ± 1.46024. The conventional 
periodic second-order scheme produces a characteristic polynomial A 2 (4 A 2 + 3) 2 = 0 for which the 
maximum eigenvalue is ^3/4* . These eigenvalues are distinctly different; the effective CFL will be 

smaller for the cyc23 than for the conventional second-order scheme. 

1 

For simplicity, the diagonal form of the cyc35 scheme (rl=160/7) will be used to illustrate the 
periodic form of the operator. The subelement As takes the form 


24 


45 

?! 

24 

44 

21 

-1/7 

-4/7 

19 

84 

0 

1/4 

1/3 

-1/8 

1/12 

-2/3 

0 

2/3 

-1/12 

1/8 

-1/3 

-1/4 

0 

11 

\i 

28 

19 

84 

4/7 

1/7 

44 

21 


Three subelements are used, and the resulting scheme is reduced by one row and column, which 
yields a matrix A p of the form 


0 

22 

21 

- 1/14 

- 2/7 

12 . 

168 

0 

0 

0 

19 

168 

2/7 

1/14 

22 

21 

ii 

24 

0 

1/4 

1/3 

- 1/8 

0 

0 

0 

0 

0 

0 

0 

1/12 

- 2/3 

0 

2/3 

- 1/12 

0 

0 

0 

0 

0 

0 

0 

1/8 

- 1/3 

- 1/4 

0 

11 

24 

0 

0 

0 

0 

0 . 

0 

0 

19 

168 

2/7 

1/14 

_22 

21 

0 

% 

- 1/14 

- 2/7 

19 

168 

0 

0 

0 

0 

0 

0 

0 

_li 

24 

0 

1/4 

1/3 

- 1/8 

0 

0 

0 

0 

0 

0 

0 

1/12 

- 2/3 

0 

2/3 

- 1/12 

0 

0 

0 

0 

0 

0 

0 

1/8 

- 1/3 

- 1/4 

0 

11 

24 

0 

0 

0 

- 12 . 

168 

0 

0 

0 


2/7 

1/14 

_22 

21 

0 

£ 

- 1/14 

- 2/7 

- 1/8 

0 

0 

0 

0 

0 

0 

0 

11 

24 

0 

1/4 

1/3 

- 1/12 

0 

0 

0 

0 

0 

0 

0 

1/12 

- 2/3 

0 

2/3 

n 

24 

0 

0 

0 

0 

0 

0 

0 

1/8 

- 1/3 

- 1/4 

0 


The matrix A p = P~ l Q p , where P is the diagonal matrix characterized by [2, 32/7, 12/7, 32/7, 
2, 32/7, 12/7, 32/7, 2, 32/7, 12/7, 32/7] on the main diagonal and Q p is entirely skew symmetric. The 
semidiscrete energy of the system defined by U T P p U will be unchanged for all time with this discretiza- 
tion and is, therefore, stable when advanced with a stable time-advancement scheme. The character- 
istic polynomial for the matrix A p is x 2 ( x 2 + 2) (9408 x 8 + 23545 x 6 + 18219 x 4 + 4302 x 2 + 243) = 0. 
The roots of this polynomial are strictly imaginary and are bounded by the points ± y/2i. Because the 
cyc35 scheme exhibits superconvergence properties, it can be compared with the conventional periodic 
fourth-order central difference expression defined on twelve grid points. The characteristic polynomial 
for the central difference case is x 2 (9x 2 + 16) (16 x 2 + 27) (48 x 2 + 49) (20736 x 4 + 19296 x 2 + 3721) = 
0. All roots are imaginary and are bounded by the points ±|i. Again, note that the cyc35 scheme 
has a slightly more restrictive CFL than the conventional central-difference method on twelve points. 
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IMAG 


EIGENVALUE SPECTRUM [CYCLO-DIFFERENCE] 



Figure 1. Eigenvalue spectra for second- to fifth-order cyclo-difference schemes with 61 gridpoints. 
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